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ABSTRACT 



Context. Interpretation of interferometric observations of rapidly rotating stars requires a good model of their surface effective tem- 
perature. Until now, laws of the form T^g oc g^^ have been used, but they are only valid for slowly rotating stars. 
Aims. We propose a simple model that can describe the latitudinal variations in the flux of rotating stars at any rotation rate. 
Methods. This model assumes that the energy flux is a divergence-free vector that is antiparallel to the effective gravity. 
Results. When mass distribution can be described by a Roche model, the latitudinal variations in the effective temperature only de- 
pend on a single parameter, namely the ratio of the equatorial velocity to the Keplerian velocity. We validate this model by comparing 
its predictions to those of the most realistic two-dimensional models of rotating stars issued from the ESTER code. The agreement is 
very good, as it is with the observations of two rapidly rotating stars, a Aql and a Leo. 

Conclusions. We suggest that as long as a gray atmosphere can be accepted, the inversion of data on flux distribution coming from 
interferometric observations of rotating stars uses such a model, which has just one free parameter. 
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1. Introduction 

The recent development of long-baseline optical/infrared inter- 
ferometry has allowed direct observation of gravity darkening at 
the surface of some rapidly rotat ing stars. This ph enomenon has 
been known since the work of von Zeipell (fT92 ?). who noticed 
that the radiative flux is proportional to the effective gravity geff 
in a barotropic star (i.e. when the pressure only depends on the 
density). Thus, the effective temperature Teff varies with latitude, 
following Teff oc g^^^. This is known as the von Zeipel law. 

Barotropicity, however, is a strong hypothesis, and is actu- 
ally incompatible with a radiative zo ne in soli d body rotation. 
This problem was noticed by E ddingtonI (Il925i) and has led to 
many discussions over subsequent decades (see lRieutordll2006l 
for a review). Now, the recent interferometri c observa t ions o f 
some nearby rapidly ro t ating stars b y McAlis ter et alJ (l2005h. 
lAufdenberg et al.' ('2006h. Ivan Belle et all (|2006). and 'Zh ao et alJ 
(.2009) have shown that gravity darkening is not well represented 
by von Zeipel's law, which seems to overestimate the tempera- 
ture difference between pole and equator. This conclusion is usu- 
ally expressed using a power law, T^ff oc g^^,, with an exponent 
less than 1 /4. This is traditionally explained by the existence of 
a thin conve ctive layer a t the surface of the star, in reference to 
the work of iLucyl ( Il967h . who found /3 ~ 0.08 for stars with a 
convective envelope. 

Actually, the weaker variation in the flux with latitude, com- 
pared to von Zeipel's law h as also been noticed by Lovekin et aL 
(2006), lEspinosa Lara & RieutordI (l2007l) . and lEspinosa Lara! 
(1201 Ol) in two-dimensional numerical models of rotating stars. 
This prompted us to reexamine this question using these new 
models. It leads us to a new approach that avoids the strong as- 
sumption of von Zeipel and that is presumably able to better 
represent the gravity darkening of fast rotating stars. This paper 
aims at presenting this new model. 



In Sect. |2] we explain the derivation of the new model of 
gravity darkening and the associated assumptions. In the next 
section, results are compared to the values of two-dimensional 
models and to observations. Conclusions follow. 



2. The model 

Before presenting our model, we briefly recall the origin of the 
von Zeipel law and Lucy's exponent. In a radiative region the 
flux is essentially carried by the diffusion of photons and thus 
correctly described by Fourier's law. 



(1) 



where is the radiative heat conductivity. If the star is assumed 
to be barotropic, density and temperature (and thus Xr) are only 
functions of the total potential O (gravitational plus centrifugal). 
Hence, 

F = -;^,(<D)vr(<D) = -xrmT'(<:»v<:> = xrmrmg^ff. (2) 

If we define the surface of the star as a place of constant 
given optical depth, we see that this model implies latitudinal 
variations in the flux that strictly follow those of the effective 
gravity. Using the effective temperature, we recover T^efF ^ 8 ■ 

Lucy's approach is based on a first-order development of 
the dependence between T^ff and geff in a convective enve- 
lope. Since the convective flux is almost orthogonal to isen- 
tropic surfaces, log Tgff and log geff are related on such a surface. 
Deep enough where the enttopy is almost the same everywhere, 
s(log Teff, loggeff) = '^0, whcrc 50 is the specific entropy in the 
bulk of the convection zone. Assuming that Tf^g oc g^^ at first 
order, one finds that 
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This equation is completed by an appropriate boundary condi- 
tion (see below). Presently, we need to note that the solution only 
depends on the effective potential, that is, on the mass distribu- 
tion. Such a quantity is generally not known, but we may observe 
that rapidly rotating stars are usually massive or intermediately 
massive, hence rather centrally condensed. We therefore chose 
the Roche model to represent the mass distribution. Its simplic- 
ity is also attractive. We thus set 



/ GM 



+ Q r sin c,. + Q r sin 6 cos f 



(7) 



where Cr and eg are radial unit vectors associated with the spher- 
ical coordinates. 

According to this model 



Fig. 1. Meridional cut representing the angle (in degrees) formed 
by the pressure gradient and the energy flux in a 3 star ro- 
tating at 90% of its Keplerian velocity at the equator calculated 
using the ESTER code. 



at zero rotation. [LucvI (1 1 967h computed this exponent for various 
stellar models and found that is was weakly dependent on the 
mass, radius, luminosity, or chemical abundance of the model. 
The mean value (3 ~ 0.08 was thus proposed. 

From the foregoing recap, we should remember that both ap- 
proaches are based on assumptions that require a small deviation 
from sphericity, and thus they demand slow rotation. 

2.1. Hypothesis for a new model 

It is quite clear from observations that slow rotation is not 
appropriate when dealing with observed gravity darkenings. 
We recall that Altair (a Aql), a favorite target of interfer- 
ometry d Pomiciano de Souza et al.ll2005l: iPeterson et alj|2006at 
iMon nier et al. 2007) is thought to be rotating close to 90% of 
the breakup angular velocity. Thus any modeling of its gravity 
darkening should be able to deal with rapid rotation. 

We propose to assume that the flux in the envelope of a ro- 
tating star is very close to 



F^-fir,e)g,«. 



(4) 



Here, / is some function of the position to be determined where 
r, 6, and (p are the spherical coordinates. In doing so, we hy- 
pothesize that the energy flux vector is almost antiparallel to the 
effective gravity. This is justified in a convective envelope where 
heat transport is essentially in the vertical direction, following 
buoyancy. In radiative regions, the flux is antiparallel to the tem- 
perature gradient, which is slightly off the potential gradient be- 
cause of baroclinicity. Fortunately, the angle between the two 
vectors remains very small, even in the case of extreme rotation. 
As shown by Fig. [1] a two-dimensional stellar model with ro- 
tation at 90% of the break-up velocity gives an angle between 
Vr and gsff of less than half a degree (i.e. less than 10"^ rd). 
Comparing our model to that of von Zeipel, we see that it does 
not imply that latitudinal variations of the flux are the same as 
those of the effective gravity or, equivalently, it allows the ratio 
F/geff to vary with latitude. 

In the envelope of a star, where no heat is generated, we have 



V-F = 



geff ■ V/ + /V ■ geff = 0. 



(5) 



(6) 



lim / (r, 0) = 



= '7 



(8) 



where L is the luminosity of the star, M its mass, and 77 is actually 
the constant that scales the function /. We therefore write 



f(r,0)^7]FJr,0) 
where 



LJ — Q. 



'cm 



(9) 



(10) 



is the nondimensional measure of the rotation rate where fl^ is 
the Keplerian angular velocity at the equator, and f - r/R^ the 
dimensionless radial coordinate scaled by the equatorial radius 
Rg. We see that in this model the latitude dependence of the flux 
is controlled by a single parameter a>. It is therefore as easy to 
use as the previous models of von Zeipel or Lucy, but it is likely 
to be much more robust at high rotation. 

2.2. Solution of the new model 

For the Roche model V ■ geff - 2Q^, so the dimensionless form 
of Eq. ^ can be written as 



1 



r sm 



df 



■ sin0cos0- 



de 



2F,., 



with the boundary condition 
F„(r = O,0) = 1. 



(11) 



(12) 



This linear, first-order, and homogeneous partial differential 
equation can be solved by using a change of variable based on 
the method of characteristics. Along a characteristic curve of 
( fTTT i we have 



dr 



A9 



-4^2 -r sin^ ( 



■ sin cos 2F,, 



(13) 



We can multiply the first and second terms by an arbitrary func- 
tion Gir, 0) 



G(r,0) 



dr 



1 



G{r,0)- 



A0 



r sm 
which leads to 



G{f,0) 



sin cos Mr ■ 



sin cos 



rsin^e de 



= dr. 



(14) 



(15) 
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where we have defined a new variable r that will be constant 
along a characteristic curve. 

Without loss of generaUty we choos43 



sin 6 



(16) 



This leads, via the chain rule dr = |^dr + ^dO, to a system of 
differential equations to calculate r 



2-2 3 

- (X) r COS ( 



St cos ( 



- cD^r" sin^cos^ 6. 



(17) 



I 86 sin6» 

From the first equation we obtain 

T = -Li??^ COS^ e + h{0) 

3 

with h(ff) an unknown function. Using (fTSl l in the second equa- 
tion of ( fTTl i. we find 



(18) 



dh cos^ G 
do ~ sine ' 

which can be integrated, yielding 



h{6) - cos + In tan - 
so that (fTSl l becomes 



T = -(w r cos^ + cos + In tan -. 
3 2 



(19) 



(20) 



(21) 



In fact, the curves r - const, are the streamlines of the flux F, 
hence also those of the effective gravity geff according to Eq. (|4|. 

After making the change of variable (r, ff) — > (r, 6), the orig- 
inal differential equation (fTTI) simplifies to (also from Eq. ( fT3]) ) 



- sm « cos f 



50 



= 2F 



which can be solved to get the general solution 

<A(T) 



tan^e' 



(22) 



(23) 



The function i^(t) can be calculated using the boundary con- 
dition (fT2l i. We define a new variable )?(r, 0) such that 

1? 1 o , , 

COS 1? -H In tan — = —cj r cos^ + cos -H In tan — - t, (24) 
2 3 2 

so that § = )?(r). For a given value of t, i? is the corresponding 
value of 6 when r — > 0. From Eq. ( l23T l we easily find that 



i/r(T) = tan^ J? 
and 

tan^j? 



tan^ e ■ 



(25) 



(26) 



' The function G(f, d) is in fact an integrating factor that should be 
chosen so that Eq. dlSt is an exact differential, namely that J; (|^) = 
then a solution T(r, ff) exists. This form is not unique since any 



This expression can be used to calculate the energy flux every- 
where in the star However, we may note singularities at = 
and - 7r/2, i.e. at the pole and the equator of the star Actually, 
these points are not real singularities since the direct resolution 
of Eq. ( fTTI ) at these points gives 



FJ0 = 0) = et' 



and 



(27) 



-2/3 



FJ0^nl2)^[l-(Jf^\"'^ (28) 
using the boundary condition (fT2t . 

2.3. The effective temperature profile 

To obtain the latitudinal profile of Teif at the surface of the star, 
one must proceed as follows. Once o) is fixed, the shape of the 
surface r{0) is given by the Roche potential 



GM 1 , , , 

+ -Q r sin - const. 

r 2 



Using Eq. ( fTOl i. it follows that r is the solution of 



+ -r^ sin^ ( 



(29) 



(30) 



for a given 0. Thus the value of )? at each point on the surface 
(r, 0) can be calculated using (l24l) . This equation is efficiently 
solved with an iterative method like the Newton one. Then we 
have 



~ [o-l ~ \47to-Gm} V tan0' 



1/4 
?eff 



/ L V7 1 4-2-2 

T -r + w sm^ ( 

\47TO-Rlj \r4 



2^2 sin^6'\''"/tan)? 



tanfl 



,(31) 



which shows that the temperature profile depends on a single 
parameter co. As a consequence the ratio of equatorial to polar 
effective temperature is, when using (l27l i and ( l28l l. given by 



^ eff 



■ eff 



-(1 - w^)'^'^exp 



4 



3 (2 -H w2)3 



(32) 



Finally we note that von Zeipel's law is recovered at slow 
rotation, since in this case & - and Eq. (|3TI ) simplifies into 



3. Results and comparison with ESTER models 

To validate the foregoing model of gravity darkening, 
we compare its results to those of fully two-dimensional 
models of r otating stars produced by the ESTEPQ code 
(e.g. Espinos a Lara & Rieutordl2007tlRieutord & Espinosa Laral 
l2009t: lEspinosa Laral I2010^ . These models utilize the OPAL 
opacities and eq uation of state (Iglesias & Rogers 199^ 
Rogers et al. 1996). They include a convective core and a fully 
radiative envelope. They are presently calculated in the limit of 



function of r will also be a solution. 



Evolution STEllaire en Rotation 
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Fig. 2. Comparison between the gravity darkening profile calculated using the model described in Sect. l2.3l (solid line), using von 
Zeipel's law (dashed line), and the models calculated using the ESTER code (crosses). Left: Star with M - 3M0 and homogeneous 
composition rotating at 70% of the Keplerian velocity at the equator (Qjt). Center. M - 3M0, Q. - O.l^k and with hydrogen 
abundance in the core Xc = Q.5X, with X the abundance in the envelope. Right: Same as center but with Q = 0.9Qt. In the bottom 
row 6 is the colatitude. 



vanishing viscosity, which requires a prescription for the differ- 
ential rotation. Here, we have chosen the surface rotation to be 
solid. The interior rotation is then self-consistently derived. 

We have computed two sets of models. The first set has 
homogeneous composition with masses between 2.5 M© and 
4 Mq and the full range of rotation velocities, from zero up to 
the breakup velocity. The second set covers the same range of 
masses and velocities, but we have set the hydrogen abundance 
in the core to be 50% of the abundance in the envelope. This 
results in a smaller core compared to the size of the star 

Figure |2] shows the comparison between our model of grav- 
ity darkening and the stellar models calculated using the ESTER 
code. It represents the profile of effective temperature as a func- 
tion of surface effective gravity (top) and colatitude (bottom). 

/ , \-l/4 

The normalized values are defined as Teff — Teff [jj^^) and 

geff - ^etf (^) ■ For comparison, we have also plotted the pro- 
file predicted by von Zeipel's law. The first two columns on the 
left show the results for a star rotating at 70% of the Keplerian 
velocity at the equatoiQ ^l^, while the third one corresponds to 
D. = Q.9D.k- The difference between the first and the second 



^ Although Keplerian velocity Q^. = ^jgelRe can be identified as the 
critical velocity, we use the former term to avoid confusion with other 
definitions of critical velocity seen in most works. A common definition 
is based on a series of models (typically Roche models) of increasing 
rotation rate, in which the critical velocity O^. is defined as the angular 
velocity of the model rotating at the break-up limit. For Roche models, 

These definitions are not equivalent, a star with 

Q. = Qi.lQ.k will have Q ^ 0.930^. 



columns is the hydrogen abundance in the core. As we observe, 
von Zeipel's law systematically overestimates the ratio between 
polar and equatorial effective temperature, while there is good 
agreement between the predictions of and the ESTER mod- 
els. The same behavior can be seen in Fig. [3] which shows the ra- 
tio between equatorial and polar effective temperature as a func- 
tion of the flattening of the star e = 1 - Rp/Re, even at high 
rotation rates. 

The good agreement between the profile of the effective 
temperature predicted by the new model and the output of the 
ESTER code is general. However, there are some slight devia- 
tions, which are more important for the models with homoge- 
neous composition. We think that this may be a consequence 
of the differential rotation of ESTER models, which contain a 
convective core rotating faster than the envelope. This leads to 
a different shape of equipotential surfaces (or isobars) compared 
to the rigidly rotating Roche model that is used in the simple 
model. This difference is less important for "evolved" models, 
because they have smaller cores. 

The gravity darkening exponent derived from a general- 
ization of von Zeipel's law T^^fi - g^j^, is commonly used in 
observational work to measure the strength of gravity darken- 
ing. In Fig.|2]we can observe that the relation between effective 
temperature and gravity is not exactly a power law. Using such 
a law is therefore not quite appropriate, but since y6-exponents 
have been derived from observations, we feel it is worth deriv- 
ing a best-fit exponent associated with ESTER models. These are 
shown in Fig.|4] This plot clearly shows that the value /3 = 0.25 
given by von Zeipel's law is appropriate only in the limit of 
slow rotation as expected. In this figure we also plotted mea- 
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4. Conclusions 




Fig. 3. Ratio between equatorial and polar effective temperature 
as a function of flattening e = 1 -Rp/Rg. Solid line: The analyti- 
cal model presented in this paper. Dashed line: Von Zeipel's law. 
Crosses: ESTER models of set 1 (Xc = X). Triangles: ESTER 
models of set 2 {X, = Q.5X). 




Observing that the energy radiated by a star is produced almost 
entirely in the stellar core and that most stars may be considered 
close to a steady state, we have noted that the energy flux is es- 
sentially a divergence-free vector field. We also observed from 
full two-dimensional models of rotating stars that the direction 
of this vector is always very close to that of the eff'ective grav- 
ity, so it only depends on the mass distribution. Following this 
picture, the physical conditions in the outer layers can only af- 
fect the flux radiated outside the star very slightly as long as a 
gray atmosphere can be assumed. Thus, we proposed that the 
energy flux in the envelope of a rotating star be approximated by 

We have shown how the non-dimensional function Fa,{r, 0) 
can be evaluated by assuming that mass distribution is repre- 
sented by the Roche model. We then demonstrated that the lati- 
tudinal variation of the effective temperature only depends on a 
single parameter u) - ^JQ?Rl/GM. Such a model is very appro- 
priate to interpreting the interferometric observations of rotating 
stars since, unlike von Zeipel's or Lucy's laws, it is valid for 
high rotation rates (up to breakup) and depends only on a single 
parameter (the yS-exponent is removed). Adjustment of the ob- 
served surface flux would thus only require variations in oj and 
/, the inclination of the rotation axis on the line of sight. 

As previously mentioned, this model fits the fully 2D mod- 
els well using a gray atmosphere with a rigidly rotating surface 
(but with interior differential rotation). This dynamical feature 
of the models may not be very realistic, so we tested a surface 
rotation Q.(9) = QgaCl - O.lcos^0) inspired by observations 
(Collier CameronI |2007|). The difference is hardly perceptible, 
therefore the proposed model of gravity darkening looks quite 
robust. Future improvement of two-dimensional models will of 
course be used to confirm this robustness. 
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Fig. 4. Gravity darkening coefficient jS as a function of flat- 
tening. The values were calculated using a least-square fit to 
the profile of effective temperature. Solid line: Model presented 
in this paper Crosses: ESTER models of set 1 {X^ - X). 
Triangles: ESTER models of set 2 (Xc - 0.5X). The coiTespond- 
ing value of /3 for von Zeipel's law is 0.25. The boxes repre- 
sent the values of /3, with their corresponding erro rs, obtained by 
inte rferometry fo r several rotating stars: Altair dMonnier et alJ 
[2007), a Cephei (iZhao et alJl2009 ). ^8 Cassiopeiae, and a Leonis 
(.Che et aL.iofll) . 



surements of the gravity darkening exponent /3 d erived by inter- 
ferometry for several rapidly rotating stars, Altair (^M onnier et alJ 
[2007), a C ephei (Zhao et al. 2009), /3 Cassiopeiae, and a Leonis 
(IChe et aril2()ll . We can see that there is good (Altair, a Leo) 
or fair (a Cep) agreement between the observed values and the 
model. The discrepancy for j6 Cas may come from its small incli- 
nation angle (~ 20 deg), so we see the star near pole-on, which 
makes determining of its flattening more difficult and makes the 
results depend more on the model used for gravity darkening and 
limb darkening. 
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